Boundary conditions for Dirac fermions on a terminated honeycomb lattice 
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We derive the boundary condition for the Dirac equation corresponding to a tight-binding model on a two- 
dimensional honeycomb lattice terminated along an arbitary direction. Zigzag boundary conditions result gener- 
ically once the boundary is not parallel to the bonds. Since a honeycomb strip with zigzag edges is gapless, 
this implies that confinement by lattice termination does not in general produce an insulating nanoribbon. We 
consider the opening of a gap in a graphene nanoribbon by a staggered potential at the edge and derive the cor- 
responding boundary condition for the Dirac equation. We analyze the edge states in a nanoribbon for arbitrary 
boundary conditions and identify a class of propagating edge states that complement the known localized edge 
states at a zigzag boundary. 

PACS numbers: 73.21.Hb, 73.22.Dj, 73.22.-f, 73.63.Bd 
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I. INTRODUCTION 

The electronic properties of graphene can be described by 
a difference equation (representing a tight-binding model on 
a honeycomb lattice) or by a differential equation (the two- 
dimensional Dirac equation) (T, ^ . The two descriptions are 
equivalent at large length scales and low energies, provided 
the Dirac equation is supplemented by boundary conditions 
consistent with the tight-binding model. These boundary con- 
ditions depend on a variety of microscopic properties, deter- 
mined by atomistic calculations |3]. 

For a general theoretical description, it is useful to know 
what boundary conditions on the Dirac equation are allowed 
by the basic physical principles of current conservation and 
(presence or absence of) time reversal symmetry — inde- 
pendently of any specific microscopic input. This problem 
was solved in Refs. 10, Ht]. The general boundary condi- 
tion depends on one mixing angle A (which vanishes if the 
boundary does not break time reversal symmetry), one three- 
dimensional unit vector n perpendicular to the normal to the 
boundary, and one three-dimensional unit vector u on the 
Bloch sphere of valley isospins. Altogether, four real parame- 
ters fix the boundary condition. 

In the present paper we investigate how the boundary condi- 
tion depends on the crystallographic orientation of the bound- 
ary. As the orientation is incremented by 30° the boundary 
configuration switches from armchair (parallel to one-third of 
the carbon-carbon bonds) to zigzag (perpendicular to another 
one-third of the bonds). The boundary conditions for the arm- 
chair and zigzag orientations are known |6]. Here we show 
that the boundary condition for intermediate orientations re- 
mains of the zigzag form, so that the armchair boundary con- 
dition is only reached for a discrete set of orientations. 

Since the zigzag boundary condition does not open up a gap 
in the excitation spectrum the implication of our result 
(not noticed in earlier studies f?*]) is that a terminated honey- 
comb lattice of arbitrary orientation is metallic rather than in- 
sulating. We present tight-binding model calculations to show 
that, indeed, the gap A oc ex-p[— f {(p)W / a] in a nanoribbon 
at crystallographic orientation (p vanishes exponentially when 
its width W becomes large compared to the lattice constant 
a, characteristic of metallic behavior. The A cx 1/W depen- 



dence characteristic of insulating behavior requires the special 
armchair orientation {(p a multiple of 60°), at which the decay 
rate f{(p) vanishes. 

Confinement by a mass term in the Dirac equation does 
produce an excitation gap regardless of the orientation of the 
boundary. We show how the infinite-mass boundary condi- 
tion of Ref. 1^ can be approached starting from the zigzag 
boundary condition, by introducing a local potential differ- 
ence on the two sublattices in the tight-binding model. Such 
a staggered potential follows from atomistic calculations 
and may well be the origin of the insulating behavior observed 
experimentally in graphene nanoribbons |9, 10]. 

The outline of this paper is as follows. In Sec.|II]we formu- 
late, following Refs. the general boundary condition of 
the Dirac equation on which our analysis is based. In Sec. Hill 
we derive from the tight-binding model the boundary condi- 
tion corresponding to an arbitrary direction of lattice termina- 
tion. In Sec.|IV]we analyze the effect of a staggered boundary 
potential on the boundary condition. In Sec. |V] we calculate 
the dispersion relation for a graphene nanoribbon with arbi- 
trary boundary conditions. We identify dispersive (= propa- 
gating) edge states which generalize the known dispersionless 
(= localized) edge states at a zigzag boundary fllll . The ex- 
ponential dependence of the gap A on the nanoribbon width 
is calculated in Sec.[yT]both analytically and numerically. We 
conclude in Sec. I VIII 



II. GENERAL BOUNDARY CONDITION 

The long-wavelength and low-energy electronic excitations 
in graphene are described by the Dirac equation 



with Hamiltonian 



H — vtq (g) (ct • p) 



(2.1) 



(2.2) 



acting on a four-component spinor wave function Here v 
is the Fermi velocity and p = ~ihV is the momentum op- 
erator. Matrices , ai are Pauli matrices in valley space and 
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sublattice space, respectively (with unit matrices tq, ctq). The 
current operator in the direction n is n • J = vtq (g) (cr • n). 

The Hamiltonian H is written in the valley isotropic rep- 
resentation of Ref. IS*]. The alternative representation H' = 
VTz IS) (cr ■ p) of Ref. |i4J is obtained by the unitary transfor- 
mation 

H' = UHU\ [/ = i(To+r,)®(7o + i(ro-T,)®a,. (2.3) 

As described in Ref. 0], the general energy-independent 
boundary condition has the form of a local linear restriction on 
the components of the spinor wave function at the boundary: 

f = M^. (2.4) 

The 4x4 matrix M has eigenvalue 1 in a two-dimensional 
subspace containing and without loss of generality we 
may assume that M has eigenvalue —1 in the orthogonal two- 
dimensional subspace. This means that M may be chosen as 
a Hermitian and unitary matrix, 

M = M\ Af2 = 1. (2.5) 

The requirement of absence of current normal to the bound- 
ary, 

{^UB ■ J|*) = 0, (2.6) 

with a unit vector normal to the boundary and pointing 
outwards, is equivalent to the requirement of anticommutation 
of the matrix M with the current operator, 

{M,nB-J} = 0. (2.7) 

That Eq. ^j} implies Eq. follows from {9\nB ■ J|*) = 
{^\M{nB ■ J)A/|*) = -(*|nB • J|^'). The converse is 
proven in App.lAl 

We are now faced with the problem of determining the most 
general 4x4 matrix M that satisfies Eqs. ( 12.51 ) and ( 12.71 ). Ref. 
[4] obtained two families of two-parameter solutions and two 
more families of three-parameter solutions. These solutions 
are subsets of the single four-parameter family of solutions 
obtained in Ref. Jl], 

M = sin A To <8) (ni • cr) + cos A (u ■ t) ig) (n2 • cr), (2.8) 

where f, ni, n2 are three-dimensional unit vectors, such that 
Til and n2 are mutually orthogonal and also orthogonal to 
ub- a proof that ( 12.81 ) is indeed the most general solution is 
given in App.|A] One can also check that the solutions of Ref. 
0] are subsets of Af = UMW. 

In this work we will restrict ourselves to boundary condi- 
tions that do not break time reversal symmetry. The time re- 
versal operator in the valley isotropic representation is 

T=-{Ty(E)ay)C, (2.9) 

with C the operator of complex conjugation. The boundary 
condition preserves time reversal symmetry if M commutes 
with T. This implies that the mixing angle A = 0, so that M 
is restricted to a three-parameter family, 

M = (v ■ t) ig {n ■ cr), n±nB. (2.10) 



III. LATTICE TERMINATION BOUNDARY 

The honeycomb lattice of a carbon monolayer is a triangu- 
lar lattice (lattice constant a) with two atoms per unit cell, re- 
ferred to as A and B atoms (see Fig.[T^). The A and B atoms 
separately form two triangular sublattices. The A atoms are 
connected only to B atoms, and vice versa. The tight-binding 
equations on the honeycomb lattice are given by 

e^Air) = t[^PB{r) + iJB{r- Ri) + V^B(r - R2)], 
etpBir) = t[i;A{r) + ^pA{r + Ri) + Va(t- + R2)]- 

Here t is the hopping energy, ipAif) and ^pB{f) are the 
electron wave functions on A and B atoms belonging to 
the same unit cell at a discrete coordinate r, while Ri = 
(a\/3/2, — a/2), R2 ~ {aV3/2, a/2) are lattice vectors as 
shown in Fig.fT^. 

Regardless of how the lattice is terminated, Eq. (13. Il l has the 
electron-hole symmetry ipB —ipB, £ — ^ ~£- For the long- 
wavelength Dirac Hamiltonian (12.2b this symmetry is trans- 
lated into the anticommutation relation 

Ha^ «) + cr^ (g) T^H = 0. (3.2) 

Electron-hole symmetry further restricts the boundary matrix 
M in Eq. (|2.10t to two classes: zigzag-like {1/ = zLz, n = z) 
and armchair-like {Vz = Uz = 0). In this section we will 
show that the zigzag-like boundary condition applies generi- 
cally to an arbitrary orientation of the lattice termination. The 
armchair-like boundary condition is only reached for special 
orientations. 



A. Characterization of the boundary 

A terminated honeycomb lattice consists of sites with three 
neighbors in the interior and sites with only one or two neigh- 
bors at the boundary. The absent neighboring sites are in- 
dicated by open circles in Fig. [T] and the dangling bonds by 
thin line segments. The tight-binding model demands that 
the wave function vanishes on the set of absent sites, so the 
first step in our analysis is the characterization of this set. We 
assume that the absent sites form a one-dimensional super- 
lattice, consisting of a supercell of N empty sites, translated 
over multiples of a superlattice vector T. Since the boundary 
superlattice is part of the honeycomb lattice, we may write 
T = nRi + mR2 with n and m non-negative integers. For 
example, in Fig. [T] we have n = 1, to = 4. Without loss 
of generality, and for later convenience, we may assume that 
m — n ~ Q (modulo 3). 

The angle ip between T and the armchair orientation (the 
X-axis in Fig.[T]) is given by 

/ 1 n — m\ n tt 

(/j = arctan — — ■ ,--<</?<-. (3.3) 

V V 3 n + ni J b b 

The armchair orientation corresponds to (p — 0, while ip — 
zLn/b corresponds to the zigzag orientation. (Because of the 
7r/3 periodicity we only need to consider \ip\ < tt/6.) 
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c) d) 

FIG. 1: (a) Honeycomb latice constructed from a unit cell (grey 
rhombus) containing two atoms (labeled A and B), translated over 
lattice vectors Ri and ^2. Panels b,c,d show three different peri- 
odic boundaries with the same period T — nRi + mi?2- Atoms on 
the boundary (connected by thick solid lines) have dangling bonds 
(thin dotted line segments) to empty neighboring sites (open circles). 
The number A*' of missing sites and A^' of dangling bonds per pe- 
riod is > n + m. Panel d shows a minimal boundary, for which 
N = N' ^n + m. 



The number N of empty sites per period T can be arbi- 
trarily large, but it cannot be smaller than n + m. Like- 
wise, the number N of dangling bonds per period cannot 
be smaller than n + m. We call the boundary minimal if 
N — N = n + m. For example, the boundary in Fig. [TJl 
is minimal {N = N =5), while the boundaries in Figs.[T]3 
and[T]: are not minimal (N ^7,N' =9 and N = 5,N' ^ 7, 
respectively). In what follows we will restrict our consider- 
ations to minimal boundaries, both for reasons of analytical 
simplicity fl?] and for physical reasons (it is natural to expect 
that the minimal boundary is energetically most favorable for 
a given orientation). 

We conclude this subsection with a property of minimal 
boundaries that we will need later on. The N empty sites 
per period can be divided into Na empty sites on sublattice 
A and Nb empty sites on sublattice B. A minimal boundary 
is constructed from n translations over Ri , each contributing 
one empty A site, and m translations over R2, each contribut- 
ing one empty B site. Hence, Na — n and Nb — m for a 
minimal boundary. 



with hk = p ■ T. While the continuous quantum number 
k E (0, 2tt) describes the propagation along the boundary, 
a second (discrete) quantum number A describes how these 
boundary modes decay away from the boundary. We select A 
by demanding that the Bloch wave ( 13.4b is also a solution of 

iPir + R3) = X4>ir). (3.5) 

The lattice vector R^ — Ri — R2 has a nonzero component 
acosip > perpendicular to T. We need |A| < 1 to 

prevent ^p{r) from diverging in the interior of the lattice. The 
decay length /decay in the direction perpendicular to T is given 
by 

-acosip 

The boundary modes satisfying Eqs. (13.4b and (13.5b are cal- 
culated in App. iBlfrom the tight-binding model. In the low- 
energy regime of interest (energies e small compared to t) 
there is an independent set of modes on each sublattice. On 
sublattice A the quantum numbers A and k are related by 

(-1-A)'"+" =exp(iA:)A" (3.7a) 

and on sublattice B they are related by 

(-1 - A)"+" = exp(ifc)A'". (3.7b) 

For a given k there are Ma roots Xp of Eq. ( I3.7ab having 
absolute value < 1, with corresponding boundary modes t/jp. 
We sort these modes according to their decay lengths from 

short to long, /decay(Ap) < /decay(Ap+l), Or |Ap| < |Ap+i|. 

The wave function on sublattice A is a superposition of these 
modes 

^fA 

(3.8) 

p=i 

with coefficients ap such that ■0'^"^) vanishes on the Na miss- 
ing A sites. Similarly there are Mb roots A^ of Eq. ( I3.7bb with 

|Ap| < 1, |Ap| < |Ap^i|. The corresponding boundary modes 
form the wave function on sublattice B, 

V'(^)=^a>;, (3.9) 
with Up such that ip'^^'> vanishes on the Nb missing B sites. 



B. Boundary modes 

The boundary breaks the two-dimensional translational in- 
variance over Ri and R2, but a one-dimensional translational 
invariance over T = nRi + mi?2 remains. The quasimo- 
mentum p along the boundary is therefore a good quantum 
number The corresponding Bloch state satisfies 

^(r + T) = exp(zA:)7/;(r), (3.4) 



C. Derivation of the boundary condition 

To derive the boundary condition for the Dirac equation it 
is sufficient to consider the boundary modes in the fc ^ 
limit. The characteristic equations ( 13.71 ) for k = each 
have a pair of solutions A± = exp(±2i7r/3) that do not 
depend on n and m. Since |A±| = 1, these modes do 
not decay as one moves away from the boundary. The cor- 
responding eigenstate exp{±iK ■ r) is a plane wave with 
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wave vector K = (4/3)7ril3/a^. One readily checks that 
this Bloch state also satisfies Eq. (|3.4| l with fc = [since 
K T ^ 27r(n - m)/3 = (modulo 27r)]. 

The wave functions ( I3.8l l and ( I3.9l l on sublattices ^ and B 
in the limit fc ^ take the form 



^^e^-f*^ '' + ^'4e-*-f*^ '" + ^ a^Vp, (3.10a) 
p=i 

A/'b-2 

p=i 

The four ampUtudes (*i, -^4) = * form the 

four-component spinor 'i' in the Dirac equation ( 12. 11 1. The re- 
maining A/a — 2 and TVs — 2 terms describe decaying boundary 
modes of the tight-binding model that are not included in the 
Dirac equation. 

We are now ready to determine what restriction on is 
imposed by the boundary condition on i/;'^"^-' and 'ip'^-^\ This 
restriction is the required boundary condition for the Dirac 
equation. In App.lBjwe calculate that, for k = 0, 



Na = n. — (n — m)/3 + 1, 
Afs = m — (m — 7i)/3 + 1, 



(3.11) 
(3.12) 



so that Ma + A/b — n + m + 2 is the total number of unknown 
amplitudes in Eqs. (13.81 1 and ( |3.9l l. These have to be chosen 
such that ijj'^^'^ and t/j^^^ vanish on Na and Nb lattice sites 
respectively. For the minimal boundary under consideration 
we have Na ~ n equations to determine A/a unknowns and 
Nb = m equations to determine As unknowns. 

Three cases can be distinguished [in each case n — m — 
(modulo 3)]: 

1. If n > m then Wa < ?i and TVs > m + 2, so 'Si = 
^'4 ~ 0, while ^'2 and are undetermined. 

2. If n < m then TVs < n and Wa > m + 2, so *2 = 
^'3 — 0, while 'i'l and ^'4 are undetermined. 



3. If n = m then A/a = n + 1 and A/s = m 
|4'i| = |*4|and|*2| = |*3|. 



1, so 



In each case the boundary condition is of the canonical form 
^ = {v> -t) ® {n- cr)* with 

1.1/ — —z,n = z if n > m (zigzag-type boundary 
condition). 

2. 1/ = z, n ^ z if n < ni (zigzag-type boundary condi- 
tion). 

3. u-z — Q,n-z = Q if n = m (armchair-type boundary 
condition). 

We conclude that the boundary condition is of zigzag-type for 
any orientation T of the boundary, unless T is parallel to the 
bonds [so that n — m and ip — (modulo 7r/3)]. 



D. Precision of the boundary condition 

At a perfect zigzag or armchair edge the four components of 
the Dirac spinor are sufficient to meet the boundary condi- 
tion. Near the boundaries with larger period and more compli- 
cated structure the wave function ( 13.101 ) also necessarily con- 
tains several boundary modes -0^, tpp that decay away from the 
boundary. The decay length S of the slowest decaying mode is 
the distance at which the boundary is indistinguishable from 
a perfect armchair or zigzag edge. At distances smaller than S 
the boundary condition breaks down. 

In the case of an armchair-like boundary (with n = to), all 
the coefficients ap and in Eqs. ( 13.10b must be nonzero to 
satisfy the boundary condition. The maximal decay length S 
is then equal to the decay length of the boundary mode tpn~i 
which has the largest |A|. It can be estimated from the char- 
acteristic equations ( 13.71 ) that S w |T|. Hence the larger the 
period of an armchair-like boundary, the larger the distance 
from the boundary at which the boundary condition breaks 
down. 

For the zigzag-like boundary the situation is different. On 
one sublattice there are more boundary modes than conditions 
imposed by the presence of the boundary and on the other 
sublattice there are less boundary modes than conditions. Let 
us assume that sublattice A has more modes than conditions 
(which happens if n < to). The quickest decaying set of 
boundary modes sufficient to satisfy the tight-binding bound- 
ary condition contains n modes t/jp with p < n. The dis- 
tance S from the boundary within which the boundary con- 
dition breaks down is then equal to the decay length of the 
slowest decaying mode 0„ in this set and is given by 



S = /decay (A„) = -O COS (/s/ 111 | A„ | . 



(3.13) 



[See Eq. diet .] 

As derived in App. IB] for the case of large periods |T| ^ 
a, the quantum number A„ satisfies the following system of 
equations; 



arg(l + A„) 



|1 + A„| 

n 



m+n 



A„| 



n + m 



arg(-A„) 



n 



n + m 



(3.14a) 
TT. (3.14b) 



The solution A„ of this equation and hence the decay length 6 
do not depend on the length |T| of the period, but only on the 
ratio n/{n + m) = (1 — V3 tan ip) /2, which is a function of 
the angle (p between T and the armchair orientation [see Eq. 
(I3.3b l. In the case n > m when sublattice B has more modes 
than conditions, the largest decay length S follows upon inter- 
changing n and to. 

As seen from Fig. |2] the resulting distance 5 within which 
the zigzag-type boundary condition breaks down is zero for 
the zigzag orientation (ip — 7r/6) and tends to infinity as the 
orientation of the boundary approaches the armchair orienta- 
tion [ip — 0). (For finite periods the divergence is cut off at 
(5 ^ |r| ^ a.) The increase of 6 near the armchair orientation 
is rather slow: For ip > 0.1 the zigzag-type boundary condi- 
tion remains precise on the scale of a few unit cells away from 
the boundary. 
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FIG. 2: Dependence on the orientation </? of the distance S from the 
boundary within which the zigzag-type boundary condition breaks 
down. The curve is calculated from formula l l3.14b valid in the limit 
\T\ 2> a of large periods. The boundary condition becomes precise 
upon approaching the zigzag orientation = n/6. 



Although the presented derivation is only valid for peri- 
odic boundaries and low energies, such that the wavelength 
is much larger than the length |T| of the boundary period, we 
argue that these conditions may be relaxed. Indeed, since the 
boundary condition is local, it cannot depend on the structure 
of the boundary far away, hence the periodicity of the bound- 
ary cannot influence the boundary condition. It can also not 
depend on the wavelength once the wavelength is larger than 
the typical size of a boundary feature (rather than the length of 
the period). Since for most boundaries both S and the scale of 
the boundary roughness are of the order of several unit cells, 
we conclude that the zigzag boundary condition is in general 
a good approximation. 



E. Density of edge states near a zigzag-like boundary 

A zigzag boundary is known to support a band of disper- 
sionless states ifTlll . which are localized within several unit 
cells near the boundary. We calculate the ID density of these 
edge states near an arbitrary zigzag-like boundary. Again as- 
suming that the sublattice A has more boundary modes than 
conditions (n < m), for each k there are AfA{k) — Na lin- 
early independent states ( 13.81 ). satisfying the boundary condi- 
tion. For fc 7^ the number of boundary modes is equal to 
Ma = n — {m — n) /3, so that for each k there are 



--AfA{k)- 



(jn — n)/3 



(3.15) 



edge states. The number of the edge states for the case when 
n > m again follows upon interchanging n and m. The den- 
sity p of edge states per unit length is given by 



P ■ 



\m — n\ 



The density of edge states is maximal p = 1 /3a for a perfect 
zigzag edge and it decreases continuously when the boundary 
orientation approaches the armchair one. Eq. (13.16b ex- 
plains the numerical data of Ref. [11], providing an analytical 
formula for the density of edge states. 



IV. STAGGERED BOUNDARY POTENTIAL 



The electron-hole symmetry ( 13.21) . which restricts the 
boundary condition to being either of zigzag-type or of 
armchair-type, is broken by an electrostatic potential. Here we 
consider, motivated by Ref. |3], the effect of a staggered po- 
tential at the zigzag boundary. We show that the effect of this 
potential is to change the boundary condition in a continuous 
way from* = ±Tz®(j2'^to * = ±T2(Ki(cr-[£xnB])*. The 
first boundary condition is of zigzag-type, while the second 
boundary condition isproduced by an infinitely large mass 
term at the boundary |i8[] . 

The staggered potential consists of a potential Va — +p, 
Vb = —p. on the ^-sites and B-sites in a total of 2N rows 
closest to the zigzag edge parallel to the y-axis (see Fig. |3]l. 
Since this potential does not mix the valleys, the boundary 
condition near a zigzag edge with staggered potential has the 
form 



^1 ^ -r^ (g) {az cos 9 + ay sin 6')^', 



(4.1) 



in accord with the general boundary condition ( 12.10b . For 9 = 
0, TT we have the zigzag boundary condition and for 6 ~ ±7r/2 
we have the infinite-mass boundary condition. 

To calculate the angle 9 we substitute Eq. ( 13.101 ) into the 
tight-binding equation ( 13.11 ) (including the staggered poten- 
tial at the left-hand side) and search for a solution in the limit 
e — 0. The boundary condition is precise for the zigzag ori- 
entation, so we may set ap = a'p = 0. It is sufficient to 
consider a single valley, so we also set ^"3 = ^4 = 0. The 
remaining nonzero components are '^le^^ '^ = ■ipA{i)e^^^ 
and 5'2e*-'^ '' = where i in the argument of iJja,b 

numbers the unit cell away from the edge and we have used 
that K points in the j/-direction. The resulting difference 
equations are 

-ptpAii) = t[ijB{i) - ^b{i - 1)], « = 1, 2, ... TV, (4.2a) 

pi^Bii) = t[^A{i) - + 1)], i = 0, 1, 2, ... TV - 1, 

(4.2b) 

V'a(0)=0. (4.2c) 

For the '$i,'i'2 components of the Dirac spinor 5* the 
boundary condition ( 14.1b is equivalent to 

^AiN)/ijBiN) = - tan(0/2). (4.3) 

Substituting the solution of Eq. (14.21 ) into Eq. ( 14.3b gives 

1 + siiih(K) sinh(K + 2NiJ./t) 



cos 9 



cosh(K) cosh(K + 2Np/t) 



(4.4) 



\T\ 3a\/ n'^ + nm + 



(3.16) 



with sinh k = p/2t. Eq. (|4.4b is exact for iV 3> 1, but it is ac- 
curate within 2% for any . The dependence of the parameter 
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9 of the boundary condition on the staggered potential strength 
/i is shown in Fig. |4]for various values of N . The boundary 
condition is closest to the infinite mass for /i/t ^ while 
the regimes ji/t <^ 1/N or /i/t S> 1 correspond to a zigzag 
boundary condition. 



2N 




i=0 1 2 3 4 



FIG. 3: Zigzag boundary with V = +fi on the j4-sites (filled dots) 
and V — — /i on the _B-sites (empty dots). The staggered poten- 
tial extends over 2N rows of atoms nearest to the zigzag edge. The 
integer i counts the number of unit cells away from the edge. 
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yU/t 

FIG. 4: Plot of the parameter 9 in the boundary condition ( i4.lt at a 
zigzag edge with the staggered potential of Fig.[3] The curves are cal- 
culated from Eq. l l4.4t . The values 6 = and 6 = 7r/2 correspond, 
respectively, to the zigzag and infinite-mass boundary conditions. 



V. DISPERSION RELATION OF A NANORIBBON 

A graphene nanoribbon is a carbon monolayer confined to 
a long and narrow strip. The energy spectrum e„(fc) of the n- 
th transverse mode is a function of the wave number k along 
the strip. This dispersion relation is nonlinear because of the 
confinement, which also may open up a gap in the spectrum 
around zero energy. We calculate the dependence of the dis- 
persion relation on the boundary conditions at the two edges 
a; = and x = W of the nanoribbon (taken along the j;-axis). 



In this section we consider the most general boundary con- 
dition ( 12. lot , constrained only by time-reversal symmetry. We 
do not require that the boundary is purely a termination of the 
lattice, but allow for arbitrary local electric fields and strained 
bonds. The conclusion of Sec.|llll that the boundary condition 
is either zigzag-like or armchair-like, does not apply therefore 
to the analysis given in this section. 

The general solution of the Dirac equation ( 12. It in the 
nanoribbon has the form '^'{x, y) = ;;(x)e''^^. We impose 
the general boundary condition (I2.10t . 

^-(0, y) = K • r) ^ (m • a)^iO, y), (5.1a) 
^{W, y) = {V2 ■ t) ® (n2 • (t)^{W, y), (5.1b) 

with three-dimensional unit vectors Ui, n^, restricted by rii ■ 
X — {i — 1, 2). (There is no restriction on the i/^.) Valley 
isotropy of the Dirac Hamiltonian (12.21 ) implies that the spec- 
trum does not depend on Ui and V2 separately but only on the 
angle 7 between them. The spectrum depends, therefore, on 
three parameters: The angle 7 and the angles 9i, 62 between 
the z-axis and the vectors rii, n2- 

The Dirac equation H'^ — e^E' has two plane wave solu- 
tions ^I* cx exp{iky + iqx) for a given e and k, corresponding 
to the two (real or imaginary) transverse wave numbers q that 
solve {hvYlk"^ + q^) = e^. Each of these two plane waves 
has a twofold valley degeneracy, so there are four indepen- 
dent solutions in total. Since the wavefunction in a ribbon is a 
linear combination of these four waves, and since each of the 
Eqs. (I5.1al5.1b] i has a two-dimensional kernel, these equations 
provide four linearly independent equations to determine four 
unknowns. The condition that Eq. ( 15. It has nonzero solutions 
gives an implicit equation for the dispersion relation of the 
nanoribbon: 

cos 6*1 cos 62 (cos uj — cos^ ri) + COS Lo sin 9i sin 62 sin^ 
— sin [sin cos 7 + sinajsin(0i — 62)] = 0, (5.2) 

where cj^ = 4W^[{e/hv)'^ — fc^] and cosfi = hvk/e. 

For 9i = 62 = and 7 = tt Eq. (5.2) reproduces the 
transcendental equation of Ref. |6| for the dispersion relation 
of a zigzag ribbon. In the case 9i =62 = tt/2 of an armchair- 
like nanoribbon, Eq. (5.2) simplifies to 

cosw = cos 7. (5.3) 

This is the only case when the transverse wave function 
'^n,k{x) is independent of the longitudinal wave number k. 
In Fig.|5]we plot the dispersion relations for several different 
boundary conditions. 

The low energy modes of a nanoribbon with |e| < hv\k\ 
[see panels a-d of Fig.|5] have imaginary transverse momen- 
tum since q^ = {e/Hv)"^ — k^ < 0. If \q\ becomes larger 
than the ribbon width W, the corresponding wave function 
becomes localized at the edges of the nanoribbon and decays 
in the bulk. The dispersion relation (15. 2t for such an edge state 
simplifies to e — hv\k \ sin 9i for the state localized near x — 
and e — -~hv\k\ sin 6*2 for the state localized near x = W. 
These dispersive edge states with velocity v sin 9 generalize 
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7=7r, 0,=Q, 0.^0 7=7r, e,=TT/5, 6l,=-7r/10 




7=0, 6i,=7r/2, 6',=7r/2 7=37r/4, 6',=7r/2, 6l2=7r/2 




-10 10 -10 10 



k [1/W] k [1/W] 

FIG. 5: Dispersion relation of nanoribbons with different boundary 
conditions. Tire large- wave number asymptotes |£| = hv\k\ of bulk 
states are shown by dashed lines. Modes that do not approach these 
asymptotes are edge states with dispersion |e| = hv\ksin6i\. The 
zigzag ribbon with 7 = tt and Oi = 62 — (a) exhibits disper- 
sionless edge states at zero energy |11]. If 6*1 or ^2 are nonzero (b, 
c) the edge states acquire linear dispersion and if sin 9i sin 82 > 
(c) a band gap opens. If 7 is unequal to or tt (d) the valleys are 
mixed which makes all the level crossings avoided and opens a band 
gap. Armchair-like ribbons with 61 — 62 — tv /2 (e, f) are the only 
ribbons having no edge states. 



the known flT'] dispersionless edge states at a zigzag bound- 
ary (with sin 9 — 0). 

Inspection of the dispersion relation i5.2\ gives the follow- 
ing condition for the presence of a gap in the spectrum of 
the Dirac equation with arbitrary boundary condition: Either 
the valleys should be mixed (7 ^ 0, tt) or the edge states 
at opposite boundaries should have energies of opposite sign 
(sin 9i sin ^2 > for 7 = tt or sin 9i sin ^2 < for 7 = 0). 

As an example, we calculate the band gap for the staggered 
potential boundary condition of Sec. |IV] We assume that the 



opposite zigzag edges have the same staggered potential, so 
that the boundary condition is 

^'(0, y) = (g) (ct^ cos 9 + ay sin 6l)^'(0, y), (5.4a) 
<ii(W,y) = -T^ ® {a:,coii9 + aySme)'^{W,y). (5.4b) 

The dependence of 9 on the parameters ji, N of the staggered 
potential is given by Eq. ( |4.4t . This boundary condition cor- 
responds to 7 = TT, 6'i = ^?2 = so that it has a gap for any 
nonzero 9. As shown in Fig.|6] A(^?) increases monotonically 
with 9 from the zigzag limit A(0) = to the infinite-mass 
limit A(7r/2) = Trhv/W. 
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FIG. 6: Dependence of the band gap A on the parameter 6 in the 
staggered potential boundary condition l l5.4t . 



VI. BAND GAP OF A TERMINATED HONEYCOMB 
LATTICE 

In this section we return to the case of a boundary formed 
purely by termination of the lattice. A nanoribbon with zigzag 
boundary condition has zero band gap according to the Dirac 
equation (Fig. [5^). According to the tight-binding equations 
there is a nonzero gap A, which however vanishes exponen- 
tially with increasing width W of the nanoribbon. We esti- 
mate the decay rate of A(W^) as follows. 

The low energy states in a zigzag-type nanoribbon are the 
hybridized zero energy edge states at the opposite boundaries. 
The energy e of such states may be estimated from the overlap 
between the edge states localized at the opposite edges, e — 
±{hv/W) exp(— W^/Zdecay)- In a perfect zigzag ribbon there 
are edge states with ^decay = (and e — 0), so that there is no 
band gap. For a ribbon with a more complicated edge shape 
the decay length of an edge state is limited by 5, the length 
within which the boundary condition breaks down (see Sec. 
|III]D). This length scale provides the analytical estimate of 
the band gap in a zigzag-like ribbon: 

A^^e->^/^ (6.1) 

with 5 given by Eqs. ( 13.131 ) and ( 13.141 ). 

The band gap of an armchair-like ribbon is 

A = {hv/W) arccos(cos7) (6.2) 
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[see Eq. ( 15.31 ) and panels e,f of Fig. |5|. Adding another row 
of atoms increases the nanoribbon width by one half of a unit 
cell and increases jhy K ■ = An/ 3, so the product AW 
in such a ribbon is an oscillatory function of W with a period 
of 1.5 unit cells. 

To test these analytical estimates, we have calculated 
A(Ty) numerically for various orientations and configura- 
tions of boundaries. As seen from Fig. |2l in ribbons with 
a non-armchair boundary the gap decays exponentially cx 
exp[—f{(p)W/a] as a function of W. Nanoribbons with the 
same orientation ip but different period |T| have the same de- 
cay rate /. As seen in Fig. [8] the decay rate obtained nu- 
merically agrees well with the analytical estimate f = a/S 
following from Eq. ( 16.1b (with 6 given as a function of (p in 
Fig. |2]l. The numerical results of Fig. |7] are consistent with 
earlier studies of the orientation dependence of the band gap 
in nanoribbons but the exponential decrease of the gap for 
non-armchair ribbons was not noticed in those studies. 

For completeness we show in Fig. |9] our numerical results 
for the band gap in an armchair-like nanoribbon {(p = 0). We 
see that the gap oscillates with a period of 1 .5 unit cells, in 
agreement with Eq. ( 16.21 ). 



1 




20 40 

W [a] 



FIG. 7: Dependence of the band gap A of zigzag-like nanoribbons on 
the width W. The curves in the left panel are calculated numerically 
from the tight-binding equations. The right panel shows the structure 
of the boundary, repeated periodically along both edges. 




7r/6 



FIG. 8: Dependence of the gap decay rate on the orientation ip of 
the boundary (defined in the inset of Fig. [2j. The dots are the fits 
to numerical results of the tight-binding equations, the solid curve is 
the analytical estimate ( 16. It . 
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FIG. 9: Dependence of the band gap A on the width W for an arm- 
chair ribbon (dashed line) and for a ribbon with a boundary of the 
same orientation but with a larger period (solid line). The curves are 
calculated numerically from the tight-binding equations. 



VII. CONCLUSION 

In summary, we have demonstrated that the zigzag-type 
boundary condition ^I^ = ir^ cr^'I' applies generically 
to a terminated honeycomb lattice. The boundary condition 
switches from the plus-sign to the minus-sign at the armchair 
orientation tf = Q (modulo vr/S), when the boundary is paral- 
lel to 1/3 of all the carbon-carbon bonds (see Fig.fTOli. 

The distance 6 from the edge within which the boundary 
condition breaks down is minimal (= 0) at the zigzag orien- 
tation (fi = 7r/6 (modulo 7r/3) and maximal at the armchair 
orientation. This is the length scale that governs the band gap 



A « {hv/W) exp{—W/6) in a nanoribbon of width W. We 
have tested our analytical results for A with the numerical so- 
lution of the tight-binding equations and find good agreement. 

While the lattice termination by itself can only produce 
zigzag or armchair-type boundary conditions, other types of 
boundary conditions can be reached by breaking the electron- 
hole symmetry of the tight-binding equations. We have con- 
sidered the effect of a staggered potential at a zigzag boundary 
(produced for example by edge magnetization [3]), and have 
calculated the corresponding boundary condition. It interpo- 
lates smoothly between the zigzag and infinite-mass boundary 
conditions, opening up a gap in the spectrum that depends on 
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the strength and range of the staggered potential. 

We have calculated the dispersion relation for arbitrary 
boundary conditions and found that the edge states which are 
dispersionless at a zigzag edge acquire a dispersion for more 
general boundary conditions. Such propagating edge states 
exist, for example, near a zigzag edge with staggered poten- 
tial. 

Our discovery that the zigzag boundary condition is generic 
explains the findings of several computer simulations iflllfisl 
[H] in which behavior characteristic of a zigzag edge was 
observed at non-zigzag orientations. It also implies that the 
mechanism of gap opening at a zigzag edge of Ref. ^ (pro- 
duction of a staggered potential by magnetization) applies 
generically to any (p ^ 0. This may explain why the band 
gap measurements of Ref. 1 10] produced results that did not 
depend on the crystallographic orientation of the nanoribbon. 




FIG. 10: These two graphene flakes (or quantum dots) both have the 
same zigzag-type boundary condition: '3/ — ±rz (8) a^'i'. The sign 
switches between + and — when the tangent to the boundary has an 
angle with the x-axis which is a multiple of 60° . 
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APPENDIX A: DERIVATION OF THE GENERAL 
BOUNDARY CONDITION (TE^ 

We first show that the anticommutation relation ( I2.7l i fol- 
lows from the current conservation requirement (12.6b . The 
current operator in the basis of eigenvectors of M has the 
block form 

^B-J={yt I). M=(J 0^). (Al) 

The Hermitian subblock X acts in the two-dimensional sub- 
space of eigenvectors of M with eigenvalue 1. To ensure that 
{'^\nB ■ J]"^) — for any 4* in this subspace it is necessary 
and sufficient that X — 0. The identity {ns ■ JY = 1 is 
equivalent to YY^ = 1 and Z = 0, hence {M, • J} = 0. 

We now show that the most general 4x4 matrix M that 
satisfies Eqs. (I2.5l l and ( I2.7l i has the 4-parameter form ( I2.8l l. 



Using only the Hermiticity of M, we have the 16-parameter 
representation 

3 

M= ^ (t, ®o-j)cy, (A2) 

with real coefficients . Anticommutation with the current 
operator brings this down to the 8-parameter form 

3 

M = ^n®{n,-(T), (A3) 

i=0 

where the n,'s are three-dimensional vectors orthogonal to 
ub- The absence of off-diagonal terms in A/^ requires that 
the vectors rii, n2, are multiples of a unit vector n which 
is orthogonal to tiq. The matrix A/ may now be rewritten as 

M = To (g) (no • cr) + (i> • r) (g) (n • cr). (A4) 

The equality A'P — 1 further demands + — 1, leading 
to the 4-parameter representation ( 12.8b after redefinition of the 
vectors. 



APPENDIX B: DERIVATION OF THE BOUNDARY MODES 

We derive the characteristic equation (13. 7b from the tight- 
binding equation ( 13.1b and the definitions of the boundary 
modes (13.4b and ( 13.5b . In the low energy limit e/t ^ o-/\T\ 
we may set e — ^ in Eq. (13. lb . so it splits into two decoupled 
sets of equations for the wave function on sublattices A and 
B: 

■ipB{r)+^Bir-Ri) + ipB{r-R2) = 0, (Bla) 
V'A(r) + ^JjAir + Ri) + iJAir + R2) = 0. (Bib) 

Substituting Ri by R2 + -R3 in these equations and using the 
definition ( 13.5b of A we express '4>{r + R2) through ipif)^ 

^B(r + i?2) = -(l + A)-Vi3(r-), (B2a) 
^A{r + R2) = -{l + \)i^A{r). (B2b) 

Eqs. ( 13.5b and ( IB2b together allow to find the boundary mode 
with a given value of A on the whole lattice: 

^B(r-+pi22+gH3) - A-?!-!- A)-PVs(r), (B3a) 
xPA{r + pR2 + qR:i) - - A)''VA(r), (B3b) 

with p and q arbitrary integers. Substituting ^(r + T) into Eq. 
(13. 4b from Eq. ( |B3b and using T = (n + m)i?2 + nR^ we 
arrive at the characteristic equation ( 13. 7b . 

We now find the roots of the Eq. (13.7b for a given k. It is 
sufficient to analyze the equation for sublattice A only since 
the calculation for sublattice B is the same after interchang- 
ing n and m. The analysis of Eq. ( |3.7ab simplifies in polar 
coordinates, 

|l + A|™+" = |Ar (B4) 
(m + n) arg(-l - X) ~ k - n arg(A) = 27r;, (B5) 
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with / = 0,±1,±2 . . .. The curve defined by Eq. (Ell is 
a contour on the complex plane around the point A = — 1 
which crosses points A± = — 1/2 ± j-\/3/2 (see Fig. [TTT ). 
The left-hand side of Eq. dBSb is a monotonic function of the 
position on this contour. If it increases by 2ttAI on the interval 
between two roots of the equation, then there are AZ — 1 roots 
inside this interval. For fc = both A_ and A+ are roots of 
the characteristic equation. So in this case the number A/a of 
roots lying inside the unit circle can be calculated from the 
increment of the left-hand side of Eq. ( IB5I ) between A_ and 
A+: 



^fA 



1 

2^ 



, 27r 27r 
[n + mj — + n— 
o o 



71 — m 
-l = n 1. 



(B6) 

Similarly, on sublattice B, we have (upon interchanging n and 

m). 



m — n 
Mb 1. 



(B7) 



The same method can be applied to calculate A,i. Since 
there are n — 1 roots on the contour defined by Eq. ( IB4b be- 



tween A„ and A* , the increment of the left-hand side of Eq. 
( |B5t between A* and A„ must be equal to 27r(ri — 1) « 27rn 
(for |T| ^ a), which immediately leads to Eq. ( 13.141 ) for A„. 
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FIG. 1 1 : Plot of the solutions of the characteristic equations ( IB4IIB5t 
for n = 5, m = 11, and fc = 0. The dots are the roots, the solid 
curve is the contour described by Eq. ( IB4t . and the dashed circles are 
unit circles with centers at and —1. 
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